install.packages("xtable")
install.packages("MASS")
install.packages("systemfit")
install.packages("locfit")  	# Used for local regression plots
install.packages("gplots")
install.packages("reporttools")
install.packages("epicalc")
install.packages("stargazer")
install.packages("ggplot2")
rm(list=ls())
set.seed(1)
library(foreign)
# library(optmatch)
library(xtable)
#require(nnet)
library(MASS)
library(systemfit)
library(locfit)		# Used for local regression plots
library(gplots)
library(reporttools)
library(epicalc)
library(stargazer)
library(ggplot2)
setwd("C:/Dropbox/MPAS2011 SURVEY/Analysis/201301_SMS_Paper/rep_sep/APSR Final manuscript/R")
# Functions
source("code/functions_final_NEWBLOCKS.r")
source("code/prepare_NEWBLOCKS.r")
# Implement Analyses
source("code/analysis_NEWBLOCKS.r")
ISSUET
A = table(DS$priority)/sum(table(DS$priority))
E = table(DS$priority[DS$ENG==1])/sum(table(DS$priority[DS$ENG==1]))
S = table(DS$priority[DS$SMS==1])/sum(table(DS$priority[DS$SMS==1]))
N = as.table(cbind(length(DS$priority),length(DS$priority[DS$ENG==1]),length(DS$priority[DS$SMS==1]) ))
Z=rbind(A, E)
T=rbind(A, S)
# TABLE OF MESSAGES
ISSUET = round(cbind(A, E, S),3)
ISSUET = rbind(ISSUET, N)
rownames(ISSUET) <- c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other", "Observations")
colnames(ISSUET) <- c("All subjects", "Engaged subjects", "SMS senders")
xtable(ISSUET)
(A-E)%*%(A-E)/2
(A-S)%*%(A-S)/2
def.par <- par(no.readonly = TRUE)
par(def.par)
ISSUET = round(cbind(A, E, S),3)
ISSUET = rbind(ISSUET, N)
ISSUET
table(DS$priority)/sum(table(DS$priority))
DS$priority
DS[1,]
DS$Priority
D$priority
rownames(DS)
colnames(DS)
# Basic price effect; lm version (reports 2 sided)
# Note higher "price" means more subsidization
# Note that I changed the code a little here since Y_REDU was not defiend
Y_REDU=cbind(DS$SMS, DS$SMS_B, DS$SMS_V)
for(y in 1:3){
print(coef(summary(lm(Y_REDU[,y]~DS$PRICE+as.factor(DS$BLOCK))))[2,], digits=3)
}
#top part of the table
MainT=cbind(
sapply(0:2, function(i) mean(DS$SMS[DS$PRICE==i])),
sapply(0:2, function(i) mean(DS$SMS_B[DS$PRICE==i])),
sapply(0:2, function(i) mean(DS$SMS_V[DS$PRICE==i])),
sapply(0:2, function(i) sum(DS$PRICE==i, na.rm=T)))
rownames(MainT) <- c("Full", "Subsidy", "Free")
colnames(MainT) <- c("Any SMS", "Public SMS", "Private SMS", "Obs")
xtable(MainT, digits=3)
# BOTTOM PART OF THE TABLE
# DEFINE REDUCED Y
Y_REDU=cbind(DS$SMS, DS$SMS_B, DS$SMS_V)
OUT = matrix(NA, 12,3)
for(y in 1:3){
OUT[1:3,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==1), DS$BLOCK, PM, C = (DS$PRICE!=2), addn=TRUE, round=3, onesided=3)  # T_FS
OUT[4:6,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=0), addn=TRUE, round=3, onesided=3)  # T_SO
OUT[7:9,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=1), addn=TRUE, round=3, onesided=3)  # T_FO
}
for(y in 1:3){
OUT[10:12,y]   <-ri_lm(Y_REDU[,y],  DS$PRICE,  DS$BLOCK, PM, addn=TRUE, round=3, onesided=3)   # Linear
}
# H4 test using SUR
msurH4     <- systemfit(list(MV = SMS_V~PRICE + as.factor(BLOCK), MB = SMS_B~PRICE + as.factor(BLOCK)), data=DS)
CM_H4	=  coef(summary(msurH4))
###coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"]
LH_H4 	= linearHypothesis(msurH4, "MB_PRICE  - MV_PRICE", test = "Chisq", na.rm=TRUE)
pH4_plus = round(LH_H4[[4]][2]/2,3)
pH4_plus
OUT= rbind(cbind(OUT, c(rep("",9),pH4_plus)), c("H2 test","","",""))
C1 = c("Subsidy vs. Full Price", "", "","Free vs. Subsidy","", "","Free vs. Full Price","", "", "Linear Trend", "","", "")
C2 = c(rep(c("ATE", paste("(p)"), "(N)"),3), c("Trend", "(p)", "(N)"), "")
OUT2 = cbind(C1,C2, OUT)
colnames(OUT2)<-c("Treatment", "Effect", "Any", "Public","Private", "H4 test")
print(xtable(OUT2, digits=3), include.rownames=FALSE)
#######################################################################################################
# Table 6 	Price Induced Flattening 1: Test of H 3.1
#######################################################################################################
OUT6 <- matrix(0, 6, 4)
DS$RICHBI = ifelse(DS$POORBI==0, 1,0)
DS$NACC = ifelse(DS$ACC==0, 1,0)
#diff in marginal effects
poverty.diff  <- lm(SMS ~ (POORBI*PRICE) + PRICE + POORBI + ACC + as.factor(BLOCK), data=DS)
poverty.diff2  <- lm(SMS ~ (RICHBI*PRICE) + PRICE + RICHBI + ACC + as.factor(BLOCK), data=DS)
s=summary(poverty.diff)
OUT6[2,1]     <- coef(poverty.diff)["PRICE"]
OUT6[2,2]     <- coef(poverty.diff2)["PRICE"]
OUT6[2,4]     <- coef(poverty.diff)["POORBI:PRICE"]
OUT6[3,1]     <- summary(poverty.diff)$coefficients["PRICE",2]
OUT6[3,2]     <- summary(poverty.diff2)$coefficients["PRICE",2]
OUT6[3,4]     <- pt(coef(s)[, 3], s$df[2], lower=FALSE)["POORBI:PRICE"]
#diff in marginal effects
acc.diff  <- lm(SMS ~ (ACC*PRICE) + PRICE + ACC + POORBI + as.factor(BLOCK), data=DS)
acc.diff2  <- lm(SMS ~ (NACC*PRICE) + PRICE + NACC + POORBI + as.factor(BLOCK), data=DS)
m=summary(acc.diff)
OUT6[5,1] <- coef(acc.diff)["PRICE"]
OUT6[5,2] <- coef(acc.diff2)["PRICE"]
OUT6[5,4] <- coef(acc.diff)["ACC:PRICE"]
OUT6[6,1]     <- summary(acc.diff)$coefficients["PRICE",2]
OUT6[6,2]     <- summary(acc.diff2)$coefficients["PRICE",2]
OUT6[6,4] <- pt(coef(m)[, 3], m$df[2], lower=FALSE)["ACC:PRICE"]
OUT6 <- round(OUT6, 3)
rows = c("", "Marginal Effect of Subsidy by Poverty","","", "Marginal Effect of Subsidy by Political Access", "")
OUT6 <- cbind(rows, OUT6)
OUT6[1, 1:5] <- c("", "Rich", "Poor", "Hypothesized Difference", "Above-Below")
OUT6[4, 1:5] <- c("", "Low Access", "High Access", "Hypothesized Difference", "Above-Below")
xtable(OUT6)
#######################################################################################################
# Table 7	  Price Induced Flattening 2: Test of H 3.2
#######################################################################################################
# Difference Marginalized v Non-Marginalized, no controls
DS$NONMARGINAL = ifelse(DS$MARGINAL==0, 1,0)
marg.diff  <- lm(SMS ~ (PRICE*MARGINAL) 	+ PRICE + MARGINAL + as.factor(BLOCK), data=DS)
marg.diff2 <- lm(SMS ~ (PRICE*NONMARGINAL) 	+ PRICE + NONMARGINAL + as.factor(BLOCK), data=DS)
v=summary(marg.diff)
q=summary(marg.diff2)
OUT7 <- matrix(NA, 3 ,2)
OUT7[1,1] <- coef(marg.diff2)["PRICE"]
OUT7[1,2] <- pt(coef(q)[, 3], q$df[2], lower=FALSE)["PRICE"]
OUT7[2,1] <- coef(marg.diff)["PRICE"]
OUT7[2,2] <- pt(coef(v)[, 3], v$df[2], lower=FALSE)["PRICE"]
OUT7[3,1] <- coef(marg.diff)["PRICE:MARGINAL"]
OUT7[3,2] <- pt(coef(v)[, 3], v$df[2], lower=FALSE)["PRICE:MARGINAL"]
OUT7 <- round(OUT7, 3)
OUT7 <- cbind(c("Marginal Effect of Price on Marginalized", "Marginal Effect of Price on Non-marginalized", "Difference"), OUT7)
xtable(OUT7)
SEND = DS$SMS==1
h(DS$MARGINAL[SEND], DS$PRICE[SEND])
bss=bsh(DS$MARGINAL[SEND],PM[SEND,])
hist(bss)
# changed to reflect treatment as subsidy (0 is full price, 1 is subsidy, 2 is complete subsidy)
h5.3tab = function(Y, X, name, PERM=PM[SEND,]) (c(
paste("Share of ",  name, " among high price senders=", round(mean(Y[X==0]),3), sep=""),
paste("n_0  =", round(length(Y[X==0]),3)),
paste("Share of ",  name, " among mid price senders= ",round(mean(Y[X==1]),3), sep=""),
paste("n_S  =", round(length(Y[X==1]),3)),
paste("Share of ",  name, " among low price senders= ",round(mean(Y[X==2]),3), sep=""),
paste("n_F  =", round(length(Y[X==2]),3)),
paste("Trend from high to low=",  round(h(Y, X),3)),
paste("(p (neg))=",mean(h(Y, X)>=bsh(Y, PERM))),
paste("(p (two sided))=",mean(abs(h(Y, X))<=abs(bsh(Y, PERM)))),
paste("(p (pos))=",mean(h(Y, X)<=bsh(Y, PERM)))
))
h5.3tab(DS$MARGINAL[SEND],DS$PRICE[SEND], "marginal respondents")
#h5.3tab(DS$POORBI[SEND],DS$PRICE[SEND], "poor respondents")
#h5.3tab(DS$WOMEN[SEND],DS$PRICE[SEND], "women")
#######################################################################################################
# Table 8	  Types of messages
#######################################################################################################
table(DS$PubPriv, DS$SMS)
#######################################################################################################
# Figure 5	Message Content and Pricing
#######################################################################################################
#treatment as subsidy (0 is full price, 1 is subsidy, 2 is complete subsidy)
counts <- as.matrix(table(DS$PRICE, DS$PubPriv))
counts <-counts[,-1] #remove uncoded messages
fullprice <- sum(DS$SMS[DS$PRICE==0 & DS$PubPriv!=0], na.rm=TRUE)
counts1 <- counts[1,]/fullprice
halfprice <- sum(DS$SMS[DS$PRICE==1 & DS$PubPriv!=0], na.rm=TRUE)
counts2 <- counts[2,]/halfprice
freeprice <- sum(DS$SMS[DS$PRICE==2 & DS$PubPriv!=0], na.rm=TRUE)
counts3 <- counts[3,]/freeprice
counts <- rbind(counts1,counts2,counts3)
pdf(file="../L/figures/PriceEffects2NEWBLOCKS.pdf", width=7, height=5.5)
label <- c("Full Price", "Partially Subsidized", "Free")
plot <- barplot(counts, main="Messaging by Type and Price",
ylim=c(0, .7),
xlab="Type of Message", ylab="Proportion of messages within each price group",
col=c("lightgrey", "darkgrey", "black"),
legend = label,
names.arg = c("Personal", "Local", "Constituency", "Group", "National"),
cex.names=1,
beside=TRUE)
dev.off()
OUT
rbind(cbind(OUT, c(rep("",9),pH4_plus)), c("H2 test","","",""))
cbind(OUT, c(rep("",9),pH4_plus)), c("H2 test","","","")
cbind(OUT, c(rep("",9),pH4_plus))
cbind(OUT, c(rep("",13),pH4_plus))
cbind(OUT, c(rep("",12),pH4_plus))
OUT= rbind(cbind(OUT, c(rep("",12),pH4_plus)), c("H2 test","","",""))
Y_REDU=cbind(DS$SMS, DS$SMS_B, DS$SMS_V)
OUT = matrix(NA, 12,3)
for(y in 1:3){
OUT[1:3,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==1), DS$BLOCK, PM, C = (DS$PRICE!=2), addn=TRUE, round=3, onesided=3)  # T_FS
OUT[4:6,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=0), addn=TRUE, round=3, onesided=3)  # T_SO
OUT[7:9,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=1), addn=TRUE, round=3, onesided=3)  # T_FO
}
for(y in 1:3){
OUT[10:12,y]   <-ri_lm(Y_REDU[,y],  DS$PRICE,  DS$BLOCK, PM, addn=TRUE, round=3, onesided=3)   # Linear
}
# H4 test using SUR
msurH4     <- systemfit(list(MV = SMS_V~PRICE + as.factor(BLOCK), MB = SMS_B~PRICE + as.factor(BLOCK)), data=DS)
CM_H4	=  coef(summary(msurH4))
###coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"]
LH_H4 	= linearHypothesis(msurH4, "MB_PRICE  - MV_PRICE", test = "Chisq", na.rm=TRUE)
pH4_plus = round(LH_H4[[4]][2]/2,3)
pH4_plus
OUTp= rbind(cbind(OUT, c(rep("",9),pH4_plus)), c("H2 test","","",""))
OUTp= rbind(cbind(OUT, c(rep("",12),pH4_plus)), c("H2 test","","",""))
OUT
OUTp= rbind(cbind(OUT, c(rep("",11),pH4_plus)), c("H2 test","","",""))
OUTp
C1 = c("Subsidy vs. Full Price", "", "","Free vs. Subsidy","", "","Free vs. Full Price","", "", "Linear Trend", "","", "")
C2 = c(rep(c("ATE", paste("(p)"), "(N)"),3), c("Trend", "(p)", "(N)"), "")
OUT2 = cbind(C1,C2, OUTp)
OUT2
pH4_plus
CM_H4
coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"]
OUT2
OUT2[11,4]
OUT2[11,6]
OUT2[11,6]<-coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"]
colnames(OUT2)<-c("Treatment", "Effect", "Any", "Public","Private", "H4 test")
OUT2
OUT2[11,6]<-round(coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"],3)
OUT2
OUT7
h5.3tab(DS$MARGINAL[SEND],DS$PRICE[SEND], "marginal respondents")
OUT6
eng.tbl
acc.tbl
par(mfrow=c(1,2))
par(xaxt="n")
colors=c("grey85","grey20")
enghist <- barplot(Z[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("Engaged", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
colors=c("grey55","grey20")
smshist <- barplot(T[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("SMS", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
dev.off()
T
############################################
# GHS
# Date: MAY 19 2014
###################################################################
#################DATA PREPARATION ################################
D=read.dta("data/GHS_Replication_NEWBLOCKS.dta", convert.factors=FALSE) # after converting the excel into Stata format.
# New variables
#D$BLOCK=matrix(NA, length(D$V1_NAME))
#    .b =sort(unique(D$V1_NAME))
#	for(i in 1:length(.b)){D$BLOCK[D$V1_NAME==.b[i]]<-i}
D$ENG   <-  (rank(D$ENGAGED)>max(rank(D$ENGAGED))*(1-mean(D$SMS[which(D$V6_SMS==1,)])))
## create rank vars
D$WEALTH_RANK  = 100-100*rank(D$POOR)/7553
D$DIST_RANK    = 100*rank(D$distance_capital)/7553
D$ENGAGED_RANK = 100*rank(D$ENGAGED)/7553
D$ACCESS_RANK  = 100*rank(D$ACCESS)/7553
D$MARG_RANK    = 100*rank(D$MARG)/7553
D$MIDDLE <- ifelse(D$WEALTH_RANK>=25 & D$WEALTH_RANK<=75, 1, 0)
# SUBSET DATA USED FOR ANALYSIS ON EXPERIMENTAL POPULATION ONLY
DS = D[D$V6_SMS==1,]
.pctile <-  1- mean(DS$SMS)
DS$ENG  <- (rank(DS$ENGAGED)>max(rank(DS$ENGAGED))*.pctile)
###################################################################
###################PERMUTATIONS####################################
# Create permutations matrix based on original randomization procedure
# Within block sampling procedure
# sample43
price.sample = function(.blocksize) ifelse(.blocksize<3,
return(sample(0:2, .blocksize)),
return(sample(c(0:2, sample(0:2,.blocksize-3))))
)
DS$block.size = sapply(DS$DICT_PCODE2011, function(i) length(DS$DICT_PCODE2011[DS$DICT_PCODE2011==i]))
A = table(DS$priority)/sum(table(DS$priority))
E = table(DS$priority[DS$ENG==1])/sum(table(DS$priority[DS$ENG==1]))
S = table(DS$priority[DS$SMS==1])/sum(table(DS$priority[DS$SMS==1]))
N = as.table(cbind(length(DS$priority),length(DS$priority[DS$ENG==1]),length(DS$priority[DS$SMS==1]) ))
Z=rbind(A, E)
T=rbind(A, S)
# TABLE OF MESSAGES
ISSUET = round(cbind(A, E, S),3)
ISSUET = rbind(ISSUET, N)
rownames(ISSUET) <- c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other", "Observations")
colnames(ISSUET) <- c("All subjects", "Engaged subjects", "SMS senders")
xtable(ISSUET)
(A-E)%*%(A-E)/2
(A-S)%*%(A-S)/2
def.par <- par(no.readonly = TRUE)
par(def.par)
par(mfrow=c(1,2))
par(xaxt="n")
colors=c("grey85","grey20")
enghist <- barplot(Z[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("Engaged", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
colors=c("grey55","grey20")
smshist <- barplot(T[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("SMS", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
dev.off()
par(mfrow=c(1,2))
par(xaxt="n")
colors=c("grey85","grey20")
enghist <- barplot(Z[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("Engaged", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
colors=c("grey55","grey20")
smshist <- barplot(T[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("SMS", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
pdf(file="../L/figures/PriorityDistNEWBLOCKS.pdf", width=13, height=8)
par(mfrow=c(1,2))
par(xaxt="n")
colors=c("grey85","grey20")
enghist <- barplot(Z[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("Engaged", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
colors=c("grey55","grey20")
smshist <- barplot(T[1:2,], beside = TRUE,
ylim = c(0, max(T)), yla="Group Share", border="blue")
legend("topright", c("SMS", "Population"), fill=colors, cex=.95)
lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
axis(1, at=seq(1, 39, by=3), labels = FALSE)
text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
dev.off()
par(def.par)
############################################
# GHS
# Date: MAY 19 2014
############################################
# install.packages("foreign")
# install.packages("xtable")
# install.packages("MASS")
# install.packages("systemfit")
# install.packages("locfit")  	# Used for local regression plots
# install.packages("gplots")
# install.packages("reporttools")
# install.packages("epicalc")
# install.packages("stargazer")
# install.packages("ggplot2")
rm(list=ls())
set.seed(1)
library(foreign)
# library(optmatch)
library(xtable)
#require(nnet)
library(MASS)
library(systemfit)
library(locfit)		# Used for local regression plots
library(gplots)
library(reporttools)
library(epicalc)
library(stargazer)
library(ggplot2)
setwd("C:/Dropbox/MPAS2011 SURVEY/Analysis/201301_SMS_Paper/rep_sep/APSR Final manuscript/R")
# Functions
source("code/functions_final_NEWBLOCKS.r")
# Prepare Variables and Create Subsetted Dataset
source("code/prepare_NEWBLOCKS.r")
# Implement Analyses
source("code/analysis_NEWBLOCKS.r")
source("code/print.all.tables.NEWBLOCK.r")
source("code/print.all.tables.NEWBLOCKS.r")
source("code/print.all.tables.NEWBLOCKS.r")
source("code/print.all.tables.NEWBLOCKS.r")
source("code/print.all.tables.NEWBLOCKS.r")
print(xtable(OUT2, digits=3), include.rownames=FALSE)
print(xtable(OUT2, digits=2), include.rownames=FALSE)
OUT
h5.3tab(DS$MARGINAL[SEND],DS$PRICE[SEND], "marginal respondents")
